#include "StdAfx.h"
#include "CalcDissimilarity.h"

namespace dhis
{
	//------------------------------------------------------------------------------
	REAL calcRegionObjectDissimilarityCriterion(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para,REAL& mergedPerim)
	{
		REAL hresult = 0.0;
		REAL hspectral = calcSpectralDissimilarityCriterion(obj1,obj2,img,para);

		mergedPerim = 0.0;
		hresult = hspectral;
		if (para->m_shapeWeight > 0.0)
		{
			hresult *= (1.0-para->m_shapeWeight);
			hresult += para->m_shapeWeight*(calcShapeDissimilarityCriterion(obj1,obj2,img,para,mergedPerim));
		}
		pround(hresult,para->m_scalePrecision);
		return hresult;
	}

	//------------------------------------------------------------------------------
	REAL calcSpectralDissimilarityCriterion(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		switch (para->m_spectralDC)
		{
		case sdcSpectralSDSBSMV:
			return sdcSQBSMSE(obj1,obj2,img,para);
		case sdcSpectralSQBSMSE:
			return sdcSDSBSMV(obj1,obj2,img,para);
		default:
			return sdcSDSBSMV(obj1,obj2,img,para);
		}
	}

	//------------------------------------------------------------------------------
	REAL calcShapeDissimilarityCriterion(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para,REAL& mergedPerim)
	{
		switch (para->m_shapeDC)
		{
		case sdcShapeCompactness:
			return shpdcCompactness(obj1,obj2,img,para,mergedPerim);
		case sdcShapeFormFactor:
			return shpdcFormFactor(obj1,obj2,img,para,mergedPerim);
		case sdcShapeAspectRatio:
			return shpdcAspectRatio(obj1,obj2,img,para,mergedPerim);
		case sdcShapeShapeIndex:
			return shpdcShapeIndex(obj1,obj2,img,para,mergedPerim);
		}
	}

	//------------------------------------------------------------------------------
	REAL sdcSDSBSMV(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL np1 = (REAL)obj1->NumPixel();
		REAL np2 = (REAL)obj2->NumPixel();
		REAL npm = np1 + np2;
		REAL h = 0.0;
		REAL dSumm,dSumsqm;
		REAL dStdDev1 = 0.0,dStdDev2 = 0.0,dStdDevm = 0.0;
		REAL dh1 = 0.0, dh2 = 0.0, dhm = 0.0;
		for (int band = 0; band < img->RasterBandCount(); ++band)
		{
			dStdDev1 = fabs( (obj1->GetSumsq(band) / np1) - (obj1->GetSum(band) * obj1->GetSum(band))/(np1*np1) );
			dStdDev1 = sqrt(dStdDev1);

			dStdDev2 = fabs( (obj2->GetSumsq(band) / np2) - (obj2->GetSum(band) * obj2->GetSum(band))/(np2*np2) );
			dStdDev2 = sqrt(dStdDev2);

			dSumm = obj1->GetSum(band) + obj2->GetSum(band);
			dSumsqm = obj1->GetSumsq(band)+ obj2->GetSumsq(band);

			dStdDevm = fabs( (dSumsqm / npm) - (dSumm * dSumm)/(npm*npm) );
			dStdDevm = sqrt(dStdDevm);

			h += fabs(npm*dStdDevm - (np1*dStdDev1+np2*dStdDev2));
		}
		return h;
	}

	//------------------------------------------------------------------------------
	REAL  sdc1Norm(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL mean1,mean2,sqdiff;
		REAL sumsqdiff = 0.0;
		REAL np1 = (REAL)obj1->NumPixel();
		REAL np2 = (REAL)obj2->NumPixel();
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			mean1 = obj1->GetSum(i)/np1;
			mean2 = obj2->GetSum(i)/np2;
			sqdiff = mean1 - mean2;
			if (sqdiff < 0.0) sqdiff = -sqdiff;
			sumsqdiff += sqdiff;
		}
		return sumsqdiff;
	}

	//------------------------------------------------------------------------------
	REAL  sdc2Norm(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL mean1,mean2,sqdiff;
		REAL sumsqdiff = 0.0;
		REAL np1 = (REAL)obj1->NumPixel();
		REAL np2 = (REAL)obj2->NumPixel();
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			mean1 = obj1->GetSum(i)/np1;
			mean2 = obj2->GetSum(i)/np2;
			sqdiff = mean1 - mean2;
			sqdiff = sqdiff*sqdiff;
			sumsqdiff += sqdiff;
		}
		return sqrt(sumsqdiff);
	}

	//------------------------------------------------------------------------------
	REAL  sdcInfinityNorm(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL mean1,mean2,sqdiff;
		REAL sumsqdiff = 0.0;
		REAL np1 = (REAL)obj1->NumPixel();
		REAL np2 = (REAL)obj2->NumPixel();
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			mean1 = obj1->GetSum(i)/np1;
			mean2 = obj2->GetSum(i)/np2;
			sqdiff = mean1 - mean2;
			if (sqdiff < 0.0) sqdiff = -sqdiff;
			if (sqdiff > sumsqdiff) sumsqdiff = sqdiff;
		}
		return sumsqdiff;
	}

	//------------------------------------------------------------------------------
	REAL  sdcSAM(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL mean1,mean2;
		REAL norm1 = 0.0, norm2 = 0.0, scalar_prod = 0.0,result=0.0;
		REAL np1 = (REAL)obj1->NumPixel();
		REAL np2 = (REAL)obj2->NumPixel();
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			mean1 = obj1->GetSum(i)/np1;
			mean2 = obj2->GetSum(i)/np2;
			norm1 += mean1*mean1;
			norm2 += mean2*mean2;
			scalar_prod += mean1*mean2;
		}
		if (norm1 != 0.0 && norm2 != 0.0)
		{
			result = scalar_prod/sqrt(norm1*norm2);
		}

		result = (HALF_P1- acos(result))/HALF_P1;
		if (result == 1.0) { return std::numeric_limits<REAL>::max(); }
		REAL x = (1.0+result)/(1.0-result);
		result = 100.0*log(x);
		return result;
	}

	//------------------------------------------------------------------------------
	REAL  sdcSID(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL mean1,mean2;
		REAL norm1 = 0.0, norm2 = 0.0,result=0.0;
		REAL np1 = (REAL)obj1->NumPixel();
		REAL np2 = (REAL)obj2->NumPixel();
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			mean1 = obj1->GetSum(i)/np1;
			mean2 = obj2->GetSum(i)/np2;
			norm1 += mean1;
			norm2 += mean2;
		}

		if (norm1 != 0.0 && norm2 != 0.0)
		{
			for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
			{
				mean1 = obj1->GetSum(i)/np1;
				mean2 = obj2->GetSum(i)/np2;
				mean1 /= norm1;
				mean2 /= norm2;

				if (mean1 != 0.0 && mean2 != 0.0)
				{
					if ((mean1/mean2) < 1.0)
						result += mean2*log(mean2/mean1) - mean1*log(mean1/mean2);
					else
						result += mean1*log(mean1/mean2) - mean2*log(mean2/mean1);
				}
			}
		}
		return result;
	}

	//------------------------------------------------------------------------------
	REAL sdcSQBSMSE(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL mean1,mean2,sqdiff;
		REAL sumsqdiff = 0.0;
		REAL np1 = (REAL)obj1->NumPixel();
		REAL np2 = (REAL)obj2->NumPixel();
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			mean1 = obj1->GetSum(i)/np1;
			mean2 = obj2->GetSum(i)/np2;
			sqdiff = mean1 - mean2;
			sqdiff = sqdiff*sqdiff;
			sumsqdiff += sqdiff;
		}
		sumsqdiff = np1*np2*sumsqdiff;
		sumsqdiff = sumsqdiff/(np1+np2);
		return sqrt(sumsqdiff);
	}

	//------------------------------------------------------------------------------
	REAL  sdcSQBMMSE(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL mean1,mean2,sqdiff;
		REAL sumsqdiff = 0.0;
		REAL np1 = (REAL)obj1->NumPixel();
		REAL np2 = (REAL)obj2->NumPixel();
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			mean1 = obj1->GetSum(i)/np1;
			mean2 = obj2->GetSum(i)/np2;
			sqdiff = mean1 - mean2;
			sqdiff = sqdiff*sqdiff;
			if (sqdiff > sumsqdiff) sumsqdiff = sqdiff;
		}
		sumsqdiff = np1*np2*sumsqdiff;
		sumsqdiff = sumsqdiff/(np1+np2);

		return sqrt(sumsqdiff);
	}

	//------------------------------------------------------------------------------
	REAL  sdcNVD(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL mean1,mean2,result;
		REAL norm1 = 0.0, norm2 = 0.0, scalar_prod = 0.0;
		REAL np1 = (REAL)obj1->NumPixel();
		REAL np2 = (REAL)obj2->NumPixel();
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			mean1 = obj1->GetSum(i)/np1;
			mean2 = obj2->GetSum(i)/np2;
			norm1 += mean1*mean1;
			norm2 += mean2*mean2;
			scalar_prod += mean1*mean2;
		}

		result = 0.0;
		if (norm1 != 0.0 && norm2 != 0.0)
		{
			result = scalar_prod/sqrt(norm1*norm2);
		}
		result = (HALF_P1-acos(result))/HALF_P1;
		if (norm2 != 0.0)
			mean1 = norm1/norm2;
		else
			mean1 = std::numeric_limits<REAL>::max();
		if (norm1 != 0.0)
			mean2 = norm2/norm1;
		else
			mean2 = std::numeric_limits<REAL>::max();
		if (mean2 < mean1)
			mean1 = mean2;
		if ((norm1 == 0.0) && (norm2 == 0.0))
			mean1 = 1.0;
		result = 1.0 - mean1*result;

		if (result == 1.0) { return std::numeric_limits<REAL>::max(); }
		REAL x = (1.0+result)/(1.0-result);
		result = 100.0*log(x);
		return result;
	}

	//------------------------------------------------------------------------------
	REAL  sdcEntropy(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL dENT = 0.0;
		REAL ro1Mean, ro2Mean, regSum, regMean, regNPix;
		REAL entropy1, entropy2;
		REAL ro1NPix = (REAL) obj1->NumPixel();
		REAL ro2NPix = (REAL) obj2->NumPixel();
		regNPix = ro1NPix + ro2NPix;
		REAL sumsqdiff = 0.0, norm1 = 0.0, norm2 = 0.0, scalar_prod = 0.0, entropy = 0.0;
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			ro1Mean = obj1->GetSum(i)/ro1NPix;
			ro2Mean = obj2->GetSum(i)/ro2NPix;

			regSum = ro1NPix*ro1Mean + ro2NPix*ro2Mean;
			regMean = regSum/regNPix;

			entropy1 = obj1->GetSum(i)*log(ro1Mean) + obj2->GetSum(i)*log(ro2Mean);
			entropy2 = obj2->GetSum(i)*log(ro2Mean) + obj1->GetSum(i)*log(ro1Mean);
			entropy1 = (entropy1 + entropy2)/2.0;
			entropy1 = (entropy1 - regSum*log(regMean));
			if (img->BandMaxValue(i) != 0.0)
			{
				entropy1 /= img->BandMaxValue(i);
			}
			entropy += entropy1;
		}
		dENT = entropy < 0.0 ? 0.0 : entropy;
		return dENT;
	}

	//------------------------------------------------------------------------------
	REAL  sdcSAR(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para)
	{
		REAL ro1Mean, ro2Mean, regNPix;
		REAL sqdiff = 0.0,sumsqdiff=0.0,dSAR = 0.0;
		REAL ro1NPix = (REAL) obj1->NumPixel();
		REAL ro2NPix = (REAL) obj2->NumPixel();
		regNPix = ro1NPix + ro2NPix;
		REAL tmp = 0.0;
		for (unsigned int i = 0; i < obj1->GetRasterCount(); ++i)
		{
			ro1Mean = obj1->GetSum(i)/ro1NPix;
			ro2Mean = obj2->GetSum(i)/ro2NPix;

			sqdiff = ro1Mean - ro2Mean;
			if (sqdiff < 0.0) sqdiff = -sqdiff;
			tmp = (ro1NPix*ro1Mean + ro2NPix*ro2Mean);
			if (tmp != 0.0) { sqdiff /= tmp; }
			sqdiff *= (ro1NPix + ro2NPix);
		}

		dSAR = sumsqdiff*sqrt((ro1NPix*ro2NPix)/(regNPix));
		return dSAR;
	}
	
	//------------------------------------------------------------------------------
	REAL  shpdcCompactness(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para,REAL& mergedPerim)
	{
		REAL h = 0.0;
		int minx,maxx;
		BoundaryLineVec mbounds;
		_mergeBoundaryline(obj1->GetBounds(),obj2->GetBounds(),mbounds,minx,maxx);
		REAL perimeterm = _calculatePerimeter(mbounds);
		REAL bhm = (mbounds.rbegin()->m_row-mbounds.begin()->m_row+1);
		REAL numpixm = obj1->NumPixel() + obj2->NumPixel();
		mergedPerim = perimeterm;

		REAL smth1 = _calculateSmoothness(obj1->Perimeter(),obj1->MaxX()-obj1->MinX()+1,obj1->MaxY()-obj1->MinY()+1);
		REAL smth2 = _calculateSmoothness(obj2->Perimeter(),obj2->MaxX()-obj2->MinX()+1,obj2->MaxY()-obj2->MinY()+1);
		REAL smthm = _calculateSmoothness(perimeterm,maxx-minx+1,bhm);

		REAL compct1 = _calculateCompactness(obj1->NumPixel(),obj1->Perimeter());
		REAL compct2 = _calculateCompactness(obj2->NumPixel(),obj2->Perimeter());
		REAL compctm = _calculateCompactness(numpixm, perimeterm);

		REAL hcompact = fabs(numpixm*compctm-(obj1->NumPixel()*compct1+obj2->NumPixel()*compct2));
		REAL hsmooth = fabs(numpixm*smthm-(obj1->NumPixel()*smth1+obj2->NumPixel()*smth2));
		h = para->m_comptWeight*hcompact + (1.0-para->m_comptWeight)*hsmooth;
		return h;
	}

	//------------------------------------------------------------------------------
	REAL  shpdcFormFactor(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para,REAL& mergedPerim)
	{
		REAL h = 0.0;
		int minx,maxx;
		BoundaryLineVec mbounds;
		_mergeBoundaryline(obj1->GetBounds(),obj2->GetBounds(),mbounds,minx,maxx);
		REAL perimeterm = _calculatePerimeter(mbounds);
		REAL bhm = (mbounds.rbegin()->m_row-mbounds.begin()->m_row+1);
		REAL numpixm = obj1->NumPixel() + obj2->NumPixel();
		mergedPerim = perimeterm;

		REAL smth1 = _calculateSmoothness(obj1->Perimeter(),obj1->MaxX()-obj1->MinX()+1,obj1->MaxY()-obj1->MinY()+1);
		REAL smth2 = _calculateSmoothness(obj2->Perimeter(),obj2->MaxX()-obj2->MinX()+1,obj2->MaxY()-obj2->MinY()+1);
		REAL smthm = _calculateSmoothness(perimeterm,maxx-minx+1,bhm);

		REAL compct1 = _calculateFormFactor(obj1->NumPixel(),obj1->Perimeter());
		REAL compct2 = _calculateFormFactor(obj2->NumPixel(),obj2->Perimeter());
		REAL compctm = _calculateFormFactor(numpixm, perimeterm);

		REAL hcompact = fabs(compctm-(compct1+compct2));
		REAL hsmooth = fabs(numpixm*smthm-(obj1->NumPixel()*smth1+obj2->NumPixel()*smth2));
		h = para->m_comptWeight*hcompact + (1.0-para->m_comptWeight)*hsmooth;
		return h;
	}

	//------------------------------------------------------------------------------
	// Aspect Ratio
	REAL  shpdcAspectRatio(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para,REAL& mergedPerim)
	{
		REAL h = 0.0;
		int minx,maxx;
		BoundaryLineVec mbounds;
		_mergeBoundaryline(obj1->GetBounds(),obj2->GetBounds(),mbounds,minx,maxx);
		REAL perimeterm = _calculatePerimeter(mbounds);
		REAL bhm = (mbounds.rbegin()->m_row-mbounds.begin()->m_row+1);
		REAL numpixm = obj1->NumPixel() + obj2->NumPixel();
		mergedPerim = perimeterm;

		REAL smth1 = _calculateSmoothness(obj1->Perimeter(),obj1->MaxX()-obj1->MinX()+1,obj1->MaxY()-obj1->MinY()+1);
		REAL smth2 = _calculateSmoothness(obj2->Perimeter(),obj2->MaxX()-obj2->MinX()+1,obj2->MaxY()-obj2->MinY()+1);
		REAL smthm = _calculateSmoothness(perimeterm,maxx-minx+1,bhm);

		REAL compct1 = _calculateAspectRatio(obj1->MaxX()-obj1->MinX()+1,obj1->MaxY()-obj1->MinY()+1);
		REAL compct2 = _calculateAspectRatio(obj2->MaxX()-obj2->MinX()+1,obj2->MaxY()-obj2->MinY()+1);
		REAL compctm = _calculateAspectRatio( (maxx+minx-1), bhm);

		REAL hcompact = fabs(numpixm*compctm-(obj1->NumPixel()*compct1+obj2->NumPixel()*compct2));
		REAL hsmooth = fabs(numpixm*smthm-(obj1->NumPixel()*smth1+obj2->NumPixel()*smth2));
		h = para->m_comptWeight*hcompact + (1.0-para->m_comptWeight)*hsmooth;
		return h;
	}

	//------------------------------------------------------------------------------
	REAL  shpdcShapeIndex(const CRegionObject* obj1,const CRegionObject* obj2,const CImageInfo* img,const CParameter* para,REAL& mergedPerim)
	{
		REAL h = 0.0;
		int minx,maxx;
		BoundaryLineVec mbounds;
		_mergeBoundaryline(obj1->GetBounds(),obj2->GetBounds(),mbounds,minx,maxx);
		REAL perimeterm = _calculatePerimeter(mbounds);
		REAL bhm = (mbounds.rbegin()->m_row-mbounds.begin()->m_row+1);
		REAL numpixm = obj1->NumPixel() + obj2->NumPixel();
		mergedPerim = perimeterm;

		REAL smth1 = _calculateSmoothness(obj1->Perimeter(),obj1->MaxX()-obj1->MinX()+1,obj1->MaxY()-obj1->MinY()+1);
		REAL smth2 = _calculateSmoothness(obj2->Perimeter(),obj2->MaxX()-obj2->MinX()+1,obj2->MaxY()-obj2->MinY()+1);
		REAL smthm = _calculateSmoothness(perimeterm,maxx-minx+1,bhm);

		REAL compct1 = _calculateShapeIndex(obj1->NumPixel(),obj1->Perimeter());
		REAL compct2 = _calculateShapeIndex(obj2->NumPixel(),obj2->Perimeter());
		REAL compctm = _calculateShapeIndex(numpixm, perimeterm);

		REAL hcompact = fabs(numpixm*compctm-(obj1->NumPixel()*compct1+obj2->NumPixel()*compct2));
		REAL hsmooth = fabs(numpixm*smthm-(obj1->NumPixel()*smth1+obj2->NumPixel()*smth2));
		h = para->m_comptWeight*hcompact + (1.0-para->m_comptWeight)*hsmooth;
		return h;
	}
}
